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The general structure of the s-wave fermionic superfluid pairing order parameter is discussed for 
systems in thermal equilibrium. We demonstrate that for finite-size systems with fixed boundary 
conditions the pairing-amplitude may always be chosen as a real function in space, in contrast to 
systems underlying periodic boundary conditions, with drastical consequences for several postulated 
Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) states. Using a simple mapping, we also investigate the 
consequences of our results for antiferromagnetic equilibrium states in a repulsive Hubbard model. 
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1. INTRODUCTION 

Unconventional superfluid pairing such as Fulde- 
Ferrell-Larkin-Ovchinnikov (FFLO)-states have become 
of great interest recently, especially in the context of ul- 
tracold fermionic quantum gases As is well-known, 
BCS-, FFLO-, breached-pair-, Sarma- and normal phases 
can occur in imbalanced two-component Fermi-mixtures, 
if there is an attractive interspecies s-wave interaction, 
which may be approximated as U(x — y) = U5(x — y) 
in the low-temperature limit The particle density, 
population imbalance and temperature are the relevant 
parameters determining which phase is actually realized, 
either in model calculations or in nature 

From a theoretical point of view, the key concept is 
the fermionic s-wave superfluid order parameter, which 
is generally described by a non-vanishing pair annihila- 
tion function A(x, x') = (^(x)^(x')). Here de- 
stroys a fermion at position x in real-space and (...) is 
the thermal average at inverse temperature /3 = l/ksT. 
In the context of Bose-Einstein condensates, the diagonal 
terms A(x) = A(x,x' = x) are interpreted as a bosonic 
pair annihilation average at position x Q . 

The special feature associated with FFLO-states is a 
spatially varying order parameter, which is either con- 
sidered to be a spatially varying complex function [l[ or 
to have sign changes as a real function of x A com- 
mon Ansatz for FFLO-states in translationally invariant 
systems is A(x) = |A| cxp(iq • x), where q ^ is the 
momentum carried by a Cooper pair in the FFLO-state 
[3, HH3] • A BCS-state with s-wave symmetry in a bal- 
anced mixture would be described by q = 0. 

FFLO-states may arise in experiments with superim- 
posed optical lattices which would theoretically be 
described by Hubbard-type lattice models, as well as in 
continuous systems Q. In this paper we discuss some 
general effects, arising from the finite system size, which 
is an immediate consequence of the presence of a trap 
potential in ultracold quantum gases. We show that a 
complex Ansatz for the superfluid order parameter is in 
general unnecessary in the presence of a confining poten- 
tial, both for continuous and for discrete systems. 



Another major experimental goal, apart from the 
search for FFLO-states, is presently the demonstration of 
antiferromagnetism in systems with superimposed opti- 
cal lattices and repulsive interaction . By mapping 
our results for superfluid states in the attractive Hub- 
bard model to the repulsive-C/ model, we show that the 
antiferromagnetic order parameter is aligned parallelly 
to a space-independent vector. This insight may signifi- 
cantly reduce the numerical effort in theoretical studies 
of antiferromagnetism in ultracold quantum gases. 

Furthermore we present numerical results for superflu- 
idity obtained in the saddle-point approximation 0, 
[l3j | . With these results we demonstrate, in particular, 
the importance of the Hartree terms in a trapped sys- 
tem, which in the literature are often neglected. We also 
present results for a Hubbard model with spin-dependent 
hopping and show that sign changes in the superfluid or- 
der parameter are not generally a feature of imbalanced 
systems, as seems to be the standard belief today. 

This paper is organized as follows. First, in section 2, 
we discuss the structure of the superfluid order parameter 
in finite systems, depending on the boundary conditions. 
By a particle-hole transformation, these results are then, 
in section 3, related to antiferromagnetic states in the 
rcpulsive-C/ Hubbard model. For the attractive-?/ Hub- 
bard model, we demonstrate in section 4 the importance 
of the Hartree terms in systems with non- vanishing mag- 
netization. Excited states containing vortices are studied 
in section 5, and we present numerical results for spin- 
dependent hopping in section 6. Finally section 7 con- 
tains a summary and our conclusions. 



2. STRUCTURE OF THE ORDER PARAMETER 

In order to determine the structure of the superfluid 
order parameter, we first consider the macroscopic de- 
scription of suprafluidity in continuum systems of finite 
extension and then discuss the consequences for finite 
discrete lattice systems. 

It is standard knowledge that the macroscopic super- 
fluid velocity and the superfluid current are described 
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by the phase </>(x) = arg[A(x)] of the order parameter 
A(x) [3,[l5|. Specifically, the superfluid current is de- 
fined as j(x) = |A(x)| 2 v(x), where v(x) = — V</?(x) is 
interpreted as the- superfluid velocity. This implies that 
in FFLO states with a spatially varying complex super- 
fluid order parameter a nonvanishing current occurs. We 
will now explain why such equilibrium solutions of FFLO 
form are suppressed in systems possessing a trapping po- 
tential, due to the finite size and the nature of the effec- 
tive boundary conditions of such systems. 

While the continuity equation for the superfluid order 
parameter, discussed above, can readily be derived for 
superfluid bosonic systems, the derivation for fcrmionic 
systems is more difficult, and the result requires careful 
interpretation. For fermions, an additional term appears 
in the continuity equation, which arises from the kinetic 
term of the Hamiltonian, as follows: 



0_ 

dt 



A(x)| 2 + V • j = — A* [V x • V x ,A(x, x')] X ' 



t- h.c. 
(1) 

Since, apparently, both arguments x and x' are required 
independently to determine the right hand side of the 
continuity equation, it is clearly more convenient to con- 
sider a combined 2<i-dimensional variable y = (x,x') T 
and focus on the full correlation function A(x,x') = 
A(y). One then arrives at a continuity equation of the 
following form: 



dt 



A(y)| 2 + V.j s = lA* 



+ h.c. 



(2) 



where we have defined the 2d-dimensional "supercurrent" 
j s analogously to the BEC current, but now as a function 
of the generalized 2c?-dimcnsional coordinates y. The 
right-hand side of ([2]) is the contribution from the in- 
teraction term and has the following explicit form for a 
contact interaction: 



UA* 



■0j(x)V> t (x)?/> t (x')V' i (x) - 

V^(x')^(x')^(x)^ t (x'))+h.c. 



(3) 



This contribution from the interaction term in the Hamil- 
tonian is very small (negligible for our purposes) for var- 
ious reasons. First, in our work we assume (as is cus- 
tomary for dilute BECs) that the interaction U is weak 
[2], so that this linear contribution in U will tend to be 
numerically small. Secondly, and more importantly, for 
superfluidity the diagonal terms A(x,x' = x) can be ex- 
pected to be dominant because of the short-ranged s- 
wave interaction, and for x' = x the interaction term in 
([3]) vanishes exactly. As a result one can to a very good 
approximation simplify @ to 



0_ 

di 



|A(y)| 5 



■3s 







(4) 



It will become clear below that (|4|), which does not hold 
rigorously but is nevertheless very accurate, simplifies the 
structure of the equilibrium solutions drastically. 

We now apply (|4]) to a condensate in equilibrium, i.e., 
to a time-independent low-temperature system without 
vortices. Accordingly we assume that the phase <p(y) 
of the order parameter has no singularities. From the 
literature we know that the curl of the superfluid vclocit 
is generally quantized along a contour C in a BEC 
i.e., 



:ity 
I, 



— <z> ds • v(x) = 27m; neN 
" Jc 



(5) 



For a condensate without any vortices (n = 0), the su- 
perfluid velocity v(x) has no singularities and, hence, 
satisfies V x v(x) = on account of (0). It is now easy 
to show from Gaufi's theorem that the superfluid order 
parameter A(y) may be chosen to be real, provided that 
the boundary conditions are fixed, since then A(y) = 
and j s = if x or x' lies on the boundary: 



= 



f dF-(<pj s ) = [ dVV-^j s ) 

JdV JV 



(6) 
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IV 


A(y) 



In the derivation of © we have used @| or Q, i.e., 
V ' Js = 0. It now follows from © that j s = for all 
values of y, and we may, therefore, choose A(x,x') € K 
for all {x, x'}. For periodic boundary conditions this 
argument fails, since the quantities mentioned above do 
not vanish at the boundaries. 



Lattice systems 

The properties of the superfluid order parameter, de- 
fined for smooth coordinates x, may be mapped onto its 
discrete analogue A(i), defined on a lattice, by the use of 
wave functions describing localized states w\ (x) at lattice 
positions i (ljj : if the trapping potential is switched off, 
the states u>i(x) would be the Wannier-functions in the 
lowest band. When an optical lattice is superimposed, 
the quantum gas is effectively described by a discretized 
second-quantized Hamiltonian of Hubbard- type form 0) , 

+ u E ■ 

i 

Here t is the nearest-neighbor hopping, V > describes 
the steepness of the parabolic trap, the parameters \i a 
(with a =t, \) describe the on-site potentials, and U is 
the on-site interaction strength, which may be chosen to 
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be attractive 17 1 or repulsive [11|. The center of the 
harmonic trap is located at i = 0. 

Our arguments imply that, also in this Hubbard- type 
model, the superfluid order parameter A(i) may always 
be chosen as a real function (provided, of course, that U is 
attractive), since due to the confining potential V a there 
is no occupation for large |i|, and therefore no superfluid 
current is present at the "border" of the system. Here we 
have assumed that u>i(x) <G R, which is exact for V = 
and should therefore be quite accurate also in the general 
case. 

Iskin and Williams 0] and Chen ct al. [HI assumed 
A(i) to be real for numerical simplicity. Our arguments 
show that there is, in fact, no need for a complex su- 
perfluid order parameter, if vortices are excluded. On 
the other hand there is a multitude of literature [j], 0- 
7] using various Ansatzes with complex forms of A(i) 
for trapped quantum gases, especially in the context of 
the local density approximation (LDA). These results are 
not consistent with our arguments and should therefore 
be carefully rechecked. 

Since these analytical results concerning the general 
structure of the order parameter rely on approxima- 
tive assumptions, we emphasize that, for the case of 
the Hubbard model (O, these results may also be con- 
firmed numerically (within the unrestricted Hartree- 
Fock approximation). Starting from periodic boundary 
conditions one may obtain complex (vortex-free) solu- 
tions of the order parameter A(i) by using Ao(i) = 
A exp(iq ■ i) as the starting point for a numerical it- 
eration. After convergence of the iteration process 
is reached, the sum over the squares of the imagi- 
nary parts SJ{A.(i)} 2 may be minimized by a global 
phase rotation as follows: A(i) — > A(i)exp(i<p) Vi. 
For periodic boundary conditions one finds (after the 
minimization): [£. 3{A(i)} 2 ]/E ; |A(i)| 2 ] « 0.3. If 
one reduces the hopping amplitudes across the periodic 
boundaries to a relative strength of 20% one obtains 
Ei9{A(i)} 2 ]/E; |A(i)| 2 ] « 0.008 and if one turns the 
hopping amplitudes across the periodic boundaries off 
(fixed boundary conditions) the imaginary part reduces 
to Ei S { A (i)} 2 ]/Eil A (i)| 2 ] « iO -5 , with otherwise 
the same parameters. This is in perfect agreement with 
our previous analytical statement: Also according to our 
numerical results, complex order parameters (without 
vortices) may occur for periodic boundary conditions, 
while, for fixed boundaries, the order parameter can be 
chosen as a real function, A(i) £ R. 



3. CONNECTION TO THE REPULSIVE-t/ 
MODEL 

The Hubbard model for attractive interaction is con- 
nected to the repulsive model (for bipartite lattices) via a 
canonical (so-called special particle-hole) transformation 



0: 



s. t -> (-1)% 



(8) 



where (— l) 1 is +1 on the A-sublattice and —1 on the 
-B-sublattice. Under this canonical transformation, the 
Hubbard Hamiltonian (|7j is transformed into 

U ^ - f H - E + K, ( 9 ) 

— U h^h^ + const , 



while simultaneously the superfluid order parameter 
A(i) = (cj^C;^) is transformed into a staggered magnetic 
order parameter in the ccy-plane (and vice versa): 

(c it c u ) o (-lY(cl t c n ) = (-lY(Si, x + iSi, y ) . (10) 

Hence, the Hubbard Hamiltonian (|7|) with a repulsive 
U is transformed into a Hamiltonian with an attrac- 
tive U and a spatially varying Zccman-term. On ac- 
count of this Zeeman term, any equilibrium state of §§§ 
has a vanishing occupation for the "down" -spin species 
(mi ~ 0) for large |i| and a full occupation for the 
"up"-spin species (n^ ~ 1). As a consequence, in the 
negativc-C/ model, the superfluid flow vanishes in the 
large- 1 i| regime. Hence, the superfluid order parameter 
A(i) may globally be chosen to be real for U < on 
the basis of the arguments presented in section 2. Trans- 
forming this real order parameter for U < back to the 
repulsivc-[/ model demonstrates that, in this case, the 
xy-antifcrromagnctic order may always be chosen along 
the x-dircction, confirming the numerical results found 
in jl3| within Hartree-Fock calculations and the results 
found in [2(| within R-DMFT calculations for particle- 
imbalanced Fermi-mixtures. 

In balanced systems the symmetry of the underlying 
Hubbard model is larger (l3^. namely SU(2) instead of 
U(l), and, therefore, antiferromagnetic order is allowed 
also along the ^-direction. For a balanced system, our 
arguments imply that the order parameter has a unique 
(globally fixed) direction. For instance a spiral struc- 
ture in the xz-plane cannot occur, since, for such a state, 
one can rotate the system around the x-axis [since the 
balanced Hamiltonian is fully 5?7(2)-invariant] to ob- 
tain an incommensurate xy-antiferromagnet, which is not 
allowed on account of our arguments presented above. 
Hence, also for balanced systems, the direction of the an- 
tiferromagnetic order parameter is globally defined. This 
result depends, of course, critically on the presence of 
a trap, which determines the boundary conditions. In- 
deed, in translationally invariant systems away from half- 
filling, incommensurate states are well-known and have 
been predicted analytically long ago [21 1. 
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(a) (b) ( c ) 

FIG. 1: Magnetization (a), |A(i)| (b) and ip(i) (c) in a vortex state with n = 1 depending on i-position on a 32 x 32 lattice. 
While |A(i)| depends only on the radial position, the magnetization is maximal in the trap center, where ip is singular. In 
contrast to a vortex-free system, the magnetization profile does not break the 7r/2-rotational symmetry of the lattice. The 
system parameters have been chosen as: U = —2.5, V = 0.025, /3 = 1000, fi = (m + AM-)/2 = 0.5 and A[i = 0.4 in units of t. 



4. ROLE OF THE HARTREE TERMS 

In sections 5 and 6 (see below) we will illustrate 
the above ideas and present numerical results for the 
attractive-C/ Hubbard model in the saddle-point approx- 
imation, which takes in general the well-known form 

n. t n u -> A(i) (c it c u + c^) - A 2 (i) (11) 

and may require the following comment. In several pre- 
vious publications the Hartree terms, which cor- 
respond to the second line in Eq. (fTTj) . were neglected for 
numerical simplicity, since they were assumed to con- 
stitute an irrelevant deformation of the trapping po- 
tential Vi 2 . This assumption, however, is only cor- 
rect in spatial regions where the magnetization vanishes 
{mi = n-j- — Hj, =0). In magnetized spatial regions 
the Hartree terms create an additional effective space- 
dependent Zeeman-term, which leads to a decrease of 
the local magnetization for U < 0. As a consequence, for 
the purposes of this paper, we must include the Hartree 
terms in our calculations, since they are of the same or- 
der of magnitude as the superfluidity terms [first line in 
Eq. (fTTj)]. In fact, we find that inclusion of the Hartree 
terms is physically highly significant and leads to results 
differing from previous literature fl8| , in which they were 
neglected. Specifically, if Hartree terms are neglected one 
finds parameter regions, in which the superfluid order 
parameter has sign changes in the tangential direction, 
while the magnetization displays a rapidly oscillating be- 
havior, as visible, e.g., in Fig. 4(g) in Ref. [l8| . However, 
inclusion of the Hartree terms leads to a simpler radial os- 
cillation, similar to the behavior visible, e.g., in Fig. 3(e) 
in Ref. In [li| it is shown that, if Hartree terms are 
neglected, the structure of the order parameter strongly 
depends on the filling in the center of the trap nc- In- 



cluding the Hartree terms we find the following results 
as compared to Ref. [IH : 

1. In the low-filling region (nc < 1) we find qualita- 
tively the same results. 

2. In the high-filling regime (nc ~ 2) the magnetiza- 
tion is smeared out by the effective Zeeman term, 
leading to a significant reduction of the peaks in 
the magnetization by a factor of 2-3. 

3. In the medium-filling regime we find qualitatively 
different behavior. Instead of the rapidly oscillating 
phase we find, depending on the filling, essentially 
the same structure as is found either in the low- or 
in the high-filling regime. 

5. THERMODYNAMIC ALLY STABLE 
VORTICES 

The derivation of the statement, that the superfluid 
order parameter can globally be assumed to be a real 
function, obviously breaks down if the superfluid veloc- 
ity is singular at any point in real-space, such as happens 
in the presence of vortices. We will now present numeri- 
cal results, obtained within the saddle-point approxima- 
tion, containing such vortices. Note that these solutions 
do not correspond to the minimum of the grand poten- 
tial and have, therefore, to be interpreted as excitations, 
where the superfluid fraction of the system carries a fi- 
nite angular momentum. Vortices tend to be numerically 
stable especially in the high-filling regime, in which the 
order parameter is small in the center of the trap. Re- 
sults are shown in Figs. [TJ [2] and [3l In Fig. Q] we show 
a solution with vortex- number n = 1, and, for compar- 
ison, we present a configuration lowering the grand po- 
tential at the same parameters in Fig. [3] Note that the 
magnetization and A(i) break the 7r/2- rotational lattice 
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symmetry in the energetically lower state, while they do 
not in the vortex state. In Fig. [2] we show a state with 
vortex-number n = 3 in a balanced mixture. 

These results show that it is possible to obtain numer- 
ically stable vortex-solutions in a system with a super- 
imposed strong optical lattice. In contrast to lattice-free 
systems Ljj , the angular motion of the condensate is not 
free in lattice systems, but is instead caused by nearest- 
neighbor tunneling processes. In addition, we have been 
able to determine the dependence of the magnetization 
distribution on the vortex quantum number n. 



6. NUMERICAL RESULTS FOR 
SPIN-DEPENDENT HOPPING 



We have shown that due to the finite size and the absence 
of vortices in the ground state the supcrfluid order pa- 
rameter may always be chosen as a real function, putting 
into question common complex Ansatzes widely spread 
in the literature. Furthermore we have investigated the 
role of the (often neglected) Hartree terms and showed 
how they qualitatively influence the structure of the su- 
pcrfluid order parameter. Finally we demonstrated that 
sign changes in the order parameter are not a general fea- 
ture of particle-imbalanced systems but rather a feature 
of the grand-canonical parameter fj^-—/j,^, again clarifying 
and generalizing previous findings. 



In addition to our analytical results regarding the na- 
ture of the superfluid order parameter, we now present 
some numerical results for a superfluid trapped Fermi- 
mixture. While we have shown that a real function A(x) 
is sufficient to describe an FFLO-state, sign changes in 
A(x) are commonly understood to characterize FFLO- 
states in a trapped system. Iskin and Williams Q and 
Chen et al. [18[ find these sign changes to be typical for 
an imbalanced superfluid mixture and interpret them as 
an FFLO state. Performing a saddle-point approxima- 
tion fl2l . [l3j for a modified attractive Hubbard Hamilto- 
nian ([7]) with spin-dependent nearest-neighbor hopping 
we do not find a sign change in A(x) to be typical for 
an imbalanced mixture. Spin-dependent hopping arises, 
e.g., if both spin species have different masses [22 1 



ra-|- 



(12) 



If hopping is spin-dependent (tf ^ t±), an imbalanced 
mixture arises naturally at ^ = /ij., since the free band- 
widths of both fermion species differ, in contrast to the 
case of spin- independent hopping (t-j. = £_[_ = £). 

Results for supcrfluid mixtures with asymmetric hop- 
ping (f-j- 7^ ti) arc shown in Fig. [4j As one can see, 
the balanced mixture has a sign-change in the super- 
fluid order parameter, while the imbalanced mixture 
has virtually no sign change (min{A} ea —0.009, while 
max{A} « 0.273). We conclude, therefore, that, for gen- 
eral hopping amplitudes t a , sign changes in A(i) arise 
from a difference in the chemical potentials f^i ^ /i^ and 
not from an imbalance itself. Furthermore we find that 
the local magnetization is not generally maximal at sites 
where the superfluid order parameter has its sign change. 



7. SUMMARY AND CONCLUSIONS 
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(a) (b) 

FIG. 2: |A(i)| (a) and <p(i) (b) in a vortex state with n = 3 depending on i-position on a 32 x 32 lattice. The system is balanced 
and therefore the magnetization is not shown. The system parameters have been chosen as: U = —2.5, V — 0.025, /3 = 1000, 
/i = (fpf + /u,j.)/2 = 0.5 and A/i = in units of t. 




FIG. 3: Magnetization (a), |A(i)| (b) and <p(i) (c) in the energetically lowest state with the same parameters as in Fig. \T\ A(i) 
is a real function. The magnetization breaks the 7r/2-rotational symmetry of the lattice. 




FIG. 4: Particle density n = (n-f + n_i)/2 (blue/dark line), population imbalance m = {n\ — n_i)/2 (orange/grey line) and 
superfluid order parameter A (green/light grey line) on a 30 x 30 square lattice depending on L-position at i y = (slice 
through the trap center). The parameters are chosen as U — —2.4, V = 0.04, ft — 50, /i = (/if + A t j.)/2 = —0.5, tf = 0.5 and 
ti — 1 in both (a) and (b). In the balanced case (a) we have chosen A/i = 0.275 in order to compensate the asymmetry arising 
from the unequal hopping terms, while in (b) we have A/i = 0. The resulting particle numbers are (a) ~ rif ~ 65 and (b) 
n_i w 56, rif ~ 71. A(i) has virtually no sign change as a function of real-space in the imbalanced case (b), while the sign 
changes in the balanced case (a). 



